﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace ODESolver
{
    public static partial class OdeSolver
    {

        static double[] cb = { 0.0, 1.0 / 2.0, 2.0 / 3.0, 1.0 / 3.0, 5.0 / 6.0, 1.0 / 6.0, 1.0 };
        static int[] hlpb = { 0, 1, 2, 3, 4, 5, 6 };
        /// <summary>
        /// 0 							
        ///	1/2						
        ///	2/9	4/9					
        ///	7/36	2/9	- 1/12				
        ///	- 35/144	- 55/36	35/48	15/8			
        ///	- 1/360	- 11/36	- 1/8	1/2	1/10		
        ///	- 41/260	22/13	43/156	- 118/39	32/195	80/39	
        /// </summary>
        static double[,] ab = { 
                           //  { 0, 0, 0, 0, 0, 0, 0 }, 
                             { 1.0 / 2.0, 0, 0, 0, 0, 0, 0 },
                             { 2.0 / 9.0, 4.0 / 9.0, 0, 0, 0, 0, 0 }, 
                             { 7.0 / 36.0, 2.0 / 9.0, -1.0 / 12.0, 0, 0, 0, 0 }, 
                             { -35.0 / 144.0, -55.0 / 36.0, 35.0 / 48.0, 15.0 / 8.0, 0, 0, 0 }, 
                             { -1.0 / 360.0, -11.0 / 36.0, -1.0 / 8.0, 1.0 / 2.0, 1.0 / 10.0, 0, 0 }, 
                             { -41.0 / 260.0, 22.0 / 13.0, 43.0 / 156.0, -118.0 / 39.0, 32.0 / 195.0, 80.0 / 39.0, 0 } };
        /// <summary>
        ///    13/200	0	11/40	11/40	4/25	4/25	13/200
        /// </summary>
        static double[] bb = { 13.0 / 200.0, 0.0, 11.0 / 40.0, 11.0 / 40.0, 4.0 / 25.0, 4.0 / 25.0, 13.0 / 200.0 };

        /// <summary>
        /// Integrate on time T by Butcher
        /// </summary>
        /// <param name="f">System of Ordinary Diff Equations</param>
        /// <param name="T">end point of integration</param>
        /// <param name="Tstep">time step</param>
        /// <param name="init">vector of initial state</param>
        /// <param name="arg">parametrs of ODESys</param>
        /// <returns></returns>
        static private double[] ButcherOnTime(ODESys f, double T, double Tstep, double[] init, double[] arg)
        {
            //throw new NotImplementedException();

            int N = init.Length;
            double[] Y = new double[N];
            double t = 0;//current t
            init.CopyTo(Y, 0);

            while (t <= T)
            {
                Y = ButcherOneStep(f, t, Tstep, Y, arg);
                t += Tstep;
            }

            return Y;
        }

        static private double[,] ButcherTrajectory(ODESys f, double T, double Tstep, double[] init, double[] arg)
        {
            //throw new NotImplementedException();
            int N = init.Length;
            double tmp = T / Tstep;
            int Nsteps = (int)Math.Truncate(tmp);
            double[,] Y = new double[Nsteps, N];
            double[] Ytemp = new double[N];
            double t = 0;
            init.CopyTo(Ytemp, 0);


            for (int i = 0; i < Nsteps; i++)
            {
                for (int j = 0; j < N; j++)
                    Y[i, j] = Ytemp[j];

                Ytemp = ButcherOneStep(f, t, Tstep, Ytemp, arg);
                t += Tstep;
            }

            return Y;
        }



        static private double[] ButcherOneStep(ODESys f,double t, double dT, double[] vector, double[] arg)
        {
            int N = vector.Length;
            double[] y = new double[N];
            int[] hlp = new int[N];
            double[] u = new double[N];
            double[][] k = new double[7][];
            for (int i = 0; i < hlp.Length; i++)
                hlp[i] = i;
            vector.CopyTo(y, 0);
            vector.CopyTo(u, 0);

            for (int n = 0; n < 6; n++)
            {
                k[n] = (from i in hlp
                        select f(i, t + cb[n] * dT, u, arg)).ToArray();

                for (int i = 0; i < N; i++)
                {
                    double tmp = (from m in hlpb
                                  select k[n][i] * ab[n, m]).Sum();
                    u[i] = y[i] + dT * tmp;
                }
            }
            k[6] = (from i in hlp
                    select f(i, t + cb[6] * dT, u, arg)).ToArray();

            for (int j = 0; j < N; j++)
                y[j] = y[j] + dT * (from i in hlpb
                                    select k[i][j] * bb[i]).Sum();

            return y;
        }

        
    }

   
}
